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Abstract 

We study the existence, stability, and mobility of fundamental discrete solitons in two- and three-dimensional nonlinear Schrodinger 
lattices with a combination of cubic self-focusing and quintic self-defocusing onsite nonlinearities. Several species of stationary 
solutions are constructed, and bifurcations linking their families are investigated using parameter continuation starting from the 
anti-continuum limit, and also with the help of a variational approximation. In particular, a species of hybrid solitons, intermediate 
^ between the site- and bond-centered types of the localized states (with no counterpart in the ID model), is analyzed in 2D and 
Ch 3D lattices. We also discuss the mobility of multi-dimensional discrete solitons that can be set in motion by lending them kinetic 
energy exceeding the appropriately crafted Peierls-Nabarro barrier; however, they eventually come to a halt, due to radiation loss. 
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1. Introduction 

A large number of models relevant to various fields of 
physics is based on discrete nonlinear Schrodinger (DNLS) 
equations [JJ. A realization of the one-dimensional (ID) 
DNLS model in arrays of parallel optical waveguides was 
predicted in Ref. [2], and later demonstrated experimen- 
tally, using an array mounted on a common substrate [3]. 
Multi-channel waveguiding systems can also be created as 
photonic lattices in bulk photorefractive crystals [JJ. Dis- 
crete solitons are fundamental self-supporting modes in the 
DNLS system [1]. The mobility 06] and collisions [BJ7] of 
discrete solitons have been studied in ID systems of the 
DNLS type with the simplest self-focusing cubic (Kerr) 
nonlinearity. The DNLS equation with the cubic onsite non- 
linearity is also a relevant model for the Bosc-Einstcin con- 
densate (BEC) trapped in deep optical lattices [§]. 
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A more general discrete cubic nonlinearity appears in the 
Salerno model [9], which combines the onsite cubic terms 
and nonlinear coupling between adjacent sites. A modifi- 
cation of the Salerno model, with opposite signs in front of 
the onsite and inter-site cubic terms, makes it possible to 
study the competition between self-focusing and defocus- 
ing discrete nonlinearities. This has been done in both ID 
[ID] and 2D [IT] settings. 

Lattice models with saturable onsite nonlinear terms 
have been studied too. The first model of that type was in- 
troduced by Vinetskii and Kukhtarev in 1975 [12]. Bright 
solitons in this model were predicted in ID [13] and 2D [14] 
geometries. Lattice solitons supported by saturable self- 
defocusing nonlinearity were created in an experiment con- 
ducted in an array of optical waveguides built in a photo- 
voltaic medium [15] . Dark discrete solitons were also con- 
sidered experimentally [IB] and theoretically [IT] in the lat- 
ter model. 

The experimental observation of optical nonlinearities 
that may be fitted by a combination of self-focusing cu- 
bic and self-defocusing quintic terms [18] suggests to study 
the dynamics of solitons in the NLS equation with cubic- 
quintic (CQ) nonlinearity. A family of stable exact soliton 
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solutions to the ID continuum NLS equation of this type is 
well known [Tj5] . The possibility to build an array of parallel 
waveguides using optical materials with the CQ nonlinear- 
ity lends relevance to the consideration of the DNLS equa- 
tion with the onsite nonlinearity of the CQ type. In particu- 
lar, this DNLS equation arises as a limit case of the contin- 
uum CQ-NLS equation which includes a periodic potential 
in the form of periodic array of rectangular channels, i.e., 
the Kronig-Pcnncy lattice. Families of stable bright soli- 
tons were found in ID [20] and 2D [21] versions of the latter 
model (the latter one with a "checkerboard" 2D potential 
supports both fundamental and vortical solitons). 

The findings of a CQ-DNLS model may also be relevant 
to the case of a self-attractive BECs confined in a 2D plane 
by a "pancake" -shaped trap combined with a sufficiently 
strong 2D optical-lattice potential (although quantum ef- 
fects such as a superfluid to Mott insulator transition are 
also relevant in the latter [22]). The condensate trapped in 
each individual potential well of this configuration is de- 
scribed by the Gross-Pitacvskii equation with an extra self- 
attractive quintic term which accounts for the deviation of 
the well's shape from one-dimensionality [23j . The tunnel- 
ing of atoms between adjacent potential wells in this set- 
ting is approximated by the linear coupling between sites 
of the respective lattice. 

The simplest stationary bright solitons, of the unstag- 
gered type (without spatial oscillations in the solitons' 
tails), have been studied in the ID version of the CQ-DNLS 
model in Ref. [21] • It was demonstrated that this class 
of solitons includes infinitely many families with distinct 
symmetries. The stability of the basic families was ana- 
lyzed, and bifurcations between them were explored in a 
numerical form, and by means of a variational approxima- 
tion (VA). Dark solitons in the same model were recently 
studied [21] and, in another very recent work, staggered 
ID bright solitons as well as the mobility of unstaggered 
ones have been investigated [26] . 

The aim of the present work is to study the existence, 
stability, and mobility of bright discrete solitons in two- 
and three-dimensional (2D and 3D) NLS lattices with the 
nonlinearity of the CQ type. As suggested by the previ- 
ous works, especially Ref. [24] , the competition of the self- 
focusing cubic and sclf-dcfocusing quintic nonlincarities in 
the setting of the discrete model may readily give rise to 
multi- stability of discrete solitons, which is not possible in 
the ordinary cubic DNLS model [27 , nor in the discrete CQ 
model where both nonlinear terms are self- focusing [28| . In 
addition to that, one may expect that the CQ model shares 
many features with those including saturable nonlinearity 
|29ll4j , such as enhanced mobility of multidimensional dis- 
crete solitons (as mentioned above, mobile discrete solitons 
can be readily found in the ID CQ-DNLS equation [26]). 

The paper is organized as follows: in the next section, 
we introduce the model and outline the method used to 
construct the multi-dimensional discrete solitons. In Sec. [3] 
we focus on stability and existence regions for 2D discrete 
solitons, and the respective bifurcations. Mobility of the 2D 



solitons on the lattice is studied in Sec. [4] Section[5]reports 
extensions of these results to 3D latices. In SeclHlwe report 
analytical results obtained by means of a VA, and Sec. [7] 
concludes the paper. 

2. The model 

In dimensionlcss form, the 2D DNLS equation with the 
onsite nonlinearity of the CQ type has the following form: 

#™,m+C'A (2 Vr l ,m+2|V'n,m| 2 l/'r l ,m-|'0n,m| 4 V'n,m =0, (1) 

where ipn,m is the complex field at site {n,, m} (the ampli- 
tude of the electromagnetic field in an optical fiber, or local 
mean-field wave function in BEC), ip = dip/dt, and C > 
is the coupling constant of the lattice model. We assume an 
isotropic medium, hence the discrete Laplacian is taken as 

A (2) 1p n , m EE Tjj n+ i. m + t/'n-l.m + V'n.m+l + 4>n,m-l ~ ^n,m- 

(2) 

The CQ nonlinearity is represented by the last two terms 
in Eq. p. 

Equation |T]) conserves two dynamical invariants: norm 
(or power, in terms of optics), 
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The conserved quantities play an important role in the anal- 
ysis of the mobility of discrete solitons, see Sec. [4] below. 

Steady state solutions are sought for in the usual form, 
ipn,m = Un,m exp(— i^ti), where fi is the real frequency, and 
the real stationary lattice field u n , m satisfies the following 
discrete equation: 
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More general solutions carrying topological charge, for 
which stationary field u n _ rn is complex, fall outside of the 
scope of the present work, and will be considered elsewhere. 

In one dimension, bright-soliton solutions of Eq. ([5]) can 
be found as homoclinic orbits of the corresponding two- 
dimensional discrete map [30j . This technique was used to 
construct ID soliton solutions to the CQ-DNLS model in 
Ref. [25- Since this method is not available in higher di- 
mensions, we construct the solutions starting from the anti- 
continuum limit, C — > 0, and perform parameter continu- 
ation to C > 0. A multidimensional version of the VA can 
also be used to construct solutions for small values of C, 
see Sec. [6] below. 

In Ref. [21], two fundamental types of solutions were 
studied: site-centered and bond-centered solitons. Each 
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Fig. 1. (Color online) Solutions to Eq. Q for (fi, C) = (-0.7,0), 
which are used as seeds to find nontrivial solutions at C > (only a 
ID slice is shown, see Fig.[3]for profiles in two dimensions). Top left: 
"Tall" (blue cross markers) and "short" (red plus markers) site-cen- 
tered solutions. Top right: "Tall" and "short" bond-centered solu- 
tions. Bottom: wider extensions of the "tall" sitc-ccntcrcd solution. 
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Fig. 3. Left (from top to bottom): Contour plots of solu- 
tions of the site-centered, bond-centered, and hybrid types for 
(H,C) = (-0.7,0.1). Right: The corresponding 3D plots. 
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Fig. 2. (Color online) The continuation to C = 0.1 of the solutions 
shown in Fig. [T] 

family of solutions was further subdivided into two sub- 
families, which represent "tall" and "short" solutions for 
given parameter values. Moreover, each sub-family con- 
tains, depending on the value of C, wider solutions that 
may be built by appending extra excited sites to the soli- 
ton. The reason for the co-existence of the tall and short 
sub-families is clearly seen in the anti-continuum. If C = 0, 
Eq. (JSJ) reduces to the following algebraic equation: 
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0. 
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^n,m i ""-n^m "-71,111 

which has at most five real solutions, viz., four nontrivial 
ones, 

U n , m = ±\l 1 ± \/l + Mi (7) 

and w n ,m = (note that these are also fixed points of the 
above-mentioned discrete map in the ID case) . Obviously, 



Eq. ([7|) gives, at most, two different non-trivial amplitudes, 
that may be continued to C > 0, giving rise to tall and 
short solitons respectively. To build wider solutions, one 
has to consider multiple sites with the nonzero field. Using 
the C — solutions as seeds, we are able to generate a large 
family of solutions in the (^, C) parameter plane, as shown 
in Figs.Q]and[2] It is found that all the wide solutions tend to 
disappear through saddle-node collisions between the tall 
and short solutions as C increases, similarly to what is the 
case for the cubic DNLS problem, as discussed in Ref. [3"T] . 

Another fundamental type of solution that arises in 
higher-dimensional lattices is a hybrid between the site- 
centered and the bond-centered solutions along the two 
spatial directions, see bottom panels in Fig. [3l This type 
of hybrid solution was considered previously in the case 
of the cubic DNLS model in Ref. [32]. We only consider 
these three symmetric types of localized states, namely the 
bond-centered, site-centered, and hybrid ones (see Fig.[3J, 
together with their intermediate asymmetric counterparts 
(see Fig.[8lc) for an example). The hybrid solution admits 
other natural variations, namely any combination of the 
various types of bond-centered solutions along one axis 
and any site-centered profile along the other. Since their 
behaviors are very similar, we consider only one such type 
of solutions. 
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Fig. 4. (Color online) (a) Power (M) versus fi for C = 0.1, and respective profiles. Bottom (from left to right): Power diagrams, for (b) 
C = 0.3, (c) C = 0.4, and (d) C = 2.0, of the bond-centered and site-centered solutions. For low values of C the co-existence of multiple 
solutions at different values of fi is obvious. The "snaking" pattern gets stretched as C increases, slowly diminishing the number of solutions 
until a single solution is left. Stable and unstable solutions are represented by solid (blue) and dashed (red) curves, respectively. 



3. The existence and stability of stationary 
solutions 

Detailed existence and stability regions of all above- 
mentioned solutions are quite intricate and particularly 
hard to detect. As described for the ID case in Ref. [H] 
and mentioned above, we expect in the 2D case the ex- 
istence of a large family of solutions at low values of C, 
which gradually annihilate, through a series of bifurca- 
tions, as C — > oo (see Ref. (31] for a detailed description of 
the termination scenaria, typically through saddle-node or 
pitchfork bifurcations, for the various families of the basic 
discrete solitons as the coupling parameter is increased). 
By plotting the power M for various types of the solutions 
(site-centered, bond-centered, and hybrid, each with vari- 
ous widths) at fixed values of C against frequency /i, it is 



possible to trace the trend followed by the solutions (see 
Fig.0|. For C = 0, the exact power for each solution can 
be found. A snake like pattern extending from fx = — 1 
to fi = exists and continues for arbitrarily large pow- 
ers. This "snaking" is also displayed for different values of 
C > in Fig. Branches of the M(/z) curve with higher 
powers correspond to wider solutions. A typical progres- 
sion observed as one follows the M(/x) curve from bottom 
(low power) to top (high power) is switching between short 
and tall solutions with gradually increasing width. For 
example, the first branch, which represents short narrow 
solutions, collides with a branch of tall narrow solutions, 
which then collides with a set of short wide solutions, and 
so on. As the coupling strength increases, the power curve 
gets stretched upward. Following the stretching, the solu- 
tions gradually vanish, until there remains a single profile. 
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Fig. 5. (Color online) Top: Pitchfork bifurcations of the bond-cen- 
tered solutions and site-centered solutions for lattice coupling con- 
stant C = 0.5. Hybrid solutions are omitted here for clarity. Bottom: 
Zoom of the bifurcation scenario depicted by the rectangular region 
in the top panel. 

Similar to what was found in cubic DNLS equation in 
Rcf. |33j the bright stationary solutions in the CQ model 
also bifurcate from plane waves (near /i f=a for the CQ 
model). It is worthwhile to highlight here the increased 
level of complexity of the relevant M(/z) curves in the 
cubic-quintic model (due to the interplay of short and tall 
solution branches) in comparison to its cubic counterpart 
of Ref. [35] , which features a single change of monotonicity 
(and correspondingly of stability) between narrow and tall 
(stable) and wide and short (unstable) solutions. 

Reference |33j provides heuristic arguments for the exis- 
tence of energy thresholds for a large class of discrete sys- 
tems with dimension higher than some critical value. This 
claim was proved in Ref. [M] for DNLS models with the 
nonlinearities of the form \ip n \ 2rT+1 '4 ) n and for chains of NLS 
equations. As can be discerned in Fig. [4] such thresholds 
also exist in the case of the cubic-quintic nonlinearity. 

In Ref. [21] a stability diagram for the discrete solitons in 
the ID model was presented in the (/i, C) plane, which gave 
a clear overview of the situation. However, in the present 
situation, the M(fi) curves for various fixed values of C, 
such as those displayed in Fig. [4] provide for a better under- 
standing of relationships between different solutions. For 
example, in the (fi, C) diagram, it would appear that the 
taller solutions cease to exist at (fi, C) m (—0.6, 0.4). How- 
ever, the respective M(fj,) curve shows that narrow and 
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Fig. 6. (Color online) Bifurcations for C = 0.4 showing that all three 
fundamental modes (site-centered, bond-centered, and hybrid) are 
all connected to each other via stability exchange with asymmetric 
solutions. Two asymmetric solutions are created where the bond- 
centered solution loses stability at the bifurcation point labeled by 
'a' in the diagram. One of these asymmetric solutions is connected 
to the hybrid solution at 'b' and the other is connected to the site- 
centered solution at point 'c'. A third type of asymmetric solution 
also emanates from the bifurcation point 'c' which is connected to 
the hybrid solution at 'd'. 

wide solitons become indistinguishable at this point, and 
deciding which solution, short or tall, is annihilated be- 
comes quite arbitrary. 

A numerical linear stability analysis was performed in 
the usual way (see Ref. [24] for details) to investigate the 
stability of each of the solution branches. As one follows a 
M(ji) curve from bottom to top, the stability is typically 
swapped around each turning point, as seen in Fig. [5] How- 
ever, the stability is not switched exactly at these points, 
as this happens via asymmetric solutions (see below). 

Similar to the ID model, a pitchfork-like bifurcation oc- 
curs between the site- and bond-centered discrete solitons. 
This is more clearly seen in Fig. [S] For C = 0.5, the bond- 
centered solution loses its stability in a neighborhood of \x « 
—0.53, and asymmetric solutions are created there. There 
are multiple asymmetric solutions in this case, but only one 
curve appears in Fig. [5] since each one is just a rotation 
of the other, hence they have the same power. The bond- 
centered solution loses its stability before the site-centered 
solution regains its stability; in fact, the site-centered soli- 
ton regains the stability exactly when the asymmetric solu- 
tions collide with it. This sort of stability exchange occurs 
throughout the M(/x) curve. The top panel of Fig. [5] shows 
two such bifurcations, with a zoom of one of them shown 
in the bottom panel. 

Figure [6] shows again a site-centered solution connected 
to a bond-centered solution but also features a connection 
of the site- and bond-centered solutions via the hybrid so- 
lution. So not only are all variations (tall, short, narrow, 
etc) within each mode connected, as shown by the snake 
like power curves, but all the fundamental modes (bond- 
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centered, site-centered, hybrid) are also connected. 

We stress that the stability regions of the above men- 
tioned fundamental modes are disjoint in regions where 
they each have roughly the same power and, unlike the ID 
model, the asymmetric solutions are unstable. These fea- 
tures can be seen in Figs. [5] and [6] Note that the multi- 
stability of symmetric solutions still occurs in this case due 
to the existence of arbitrarily wide solutions at fixed values 
of C (see Fig. [4}. As a general comment, it should be noted 
that many of the features of the 2D cubic-quintic model 
(such as e.g., the existence of unstable asymmetric solu- 
tions, and their connecting the fundamental modes) can 
also be observed in the case of the saturable model [2] , al- 
though in the present case of the cubic-quintic model, the 
relevant phenomenology is even richer due to, for instance, 
the existence of multiple (i.e., tall and short) steady states. 

4. The mobility 

In one dimension, traveling solutions can be found in the 
form 

^ n = u(n-vt)e i » t , (8) 

where v is a real velocity. Substitution of this expression 
in the ID DNLS model yields the following advance-delay 
differential equation 

= -i[vu(z) + i^u(z)} + 2\u(z)\ 2 u(z) - \u{z)\ 4 u(z) 

+C [u(z + 1) + u{z - 1) - 2u(z)] , (9) 

where z = n—vt. Stationary solutions arc said to be transla- 
tionally invariant if the function u n = u(nh), where h is the 
lattice spacing, can be extended to a one-parameter fam- 
ily of continuous solutions, u(z — s), of the advance-delay 
equation ^ with v = 0. Solutions of this type have been 
found in other lattice models (see Ref. [35136] and references 
therein). Localized solutions with non-oscillatory tails in 
similar models for v ^ 0, have been found in Ref. [29 37J 
by solving a respective counterpart of Eq. ((9]). If transla- 
tionally invariant solutions exist, then the sundry modes 
(bond-centered, site-centered, etc.) are generated by the 
same continuous function u(z — s), each with a correspond- 
ing value of s. The translationally invariant solutions occur 

(i) at transparency points, which are points in the parame- 
ter space where solutions exchange their stability, and (ii) 
if the Peierls-Nabarro (PN) barrier vanishes; the barrier 
being defined as the difference in energy between the site- 
centered and bond-centered solutions. Note that (i) and 

(ii) are necessary but not sufficient conditions for the ex- 
istence of translationally invariant solutions. For higher- 
dimensional lattices, translationally invariant solutions for 
DNLS- like models have not been found yet. However, effec- 
tively mobile lattice solitons have been found in 2D models 
in regions of the parameter space where the PN barrier is 
low (enhanced mobility). This has been the case both for 
quadratic nonlinearities [35] and in the vicinity of stability 
exchanges for saturable models [14]. The resulting mobile 
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Fig. 7. (Color online) Resulting density plots, of a one-dimensional 
slice along the axis of propagation, from imprinting momentum to 
a stationary soliton by means of the "kick" defined in Eq. I10I (a) 
k < fcdepin : the solution remains pinned at its initial position, (b) 
^depin < k n < ^disperse 1 the solution becomes mobile, but eventually 
comes to a halt due to radiation loss, (c) k > ^disperse 1 the kick is 
so strong that the solution disperses. 

solutions radiate energy and eventually come to a halt. Ex- 
act solutions of the corresponding advance-delay differen- 
tial equation, if they exist, would experience no radiation 
losses and propagate indefinitely, which is why they are 
called radiationless solutions [55] . As mentioned above for 
translationally invariant solutions, radiationless solutions 
have also not been yet been found in higher-dimensional 
lattices. In fact, it is an important open question whether 
such solutions exist typically, since the single tail resonance 
appropriately made to vanish in Ref. [29^ to obtain such ex- 
ponentially localized traveling solutions in ID settings, ac- 
quires infinite multiplicity in higher dimensional settings. 
Thus, the admittedly straightforward technique of identi- 
fying regions of enhanced mobility may be the only possible 
method in higher dimensional DNLS problems. 

The goal is to "kick" the stationary solutions into mo- 
tion. From a Hamiltonian point of view, the real part of the 
solution corresponds to position and the imaginary part to 
momentum [27| . Therefore, in order to set it into motion 
one should apply a perturbation that will alter the imagi- 
nary part of the solution in an asymmetric way, and thus 
providing it with the necessary momentum to move. To this 
end, we apply a "kick" of the form 

«n,m(0) = U n , m e« k " n+k ™ m \ (10) 

where u n _ m is a stationary solution, and k n , and k m are 
real wavenumbcrs. This method has been used in numer- 
ous studies in one-dimensional settings (cf. Rcfs. |27I39| ) 
and recently in two-dimensions [14] . Bright mobile solutions 
were studied in this way in the ID CQ model in Ref. [28] and 
in greater detail (and for staggered solutions) in Ref. [26] . 
We present here results for a site-centered solution mov- 
ing along a single axis only. Therefore we set k„ ^ and 
k m = 0. Results for other solutions are similar. There are 
three qualitative scenarios that we have observed as re- 
sult of the kick: (a) the kick is below some threshold value, 
k n < ^depin, and so the corresponding energy increase is too 
low to depin the solution, (b) the kick is greater than this 
threshold value, k n > fcdepin, an d the solution is set in mo- 
tion eventually halting, or (c) the initial kick is so strong, 
k n > ^disperse , that the solution disperses. See Fig. [7] for 
examples of these three scenaria. 

We are interested in areas of parameter space that pro- 
vide good conditions for mobility for the kicked solutions. 
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Fig. 8. (Color online) Evolution of a site-centered soliton kicked along 
a diagonal. In the course of its motion, the traveling object takes on 
the (a) site-centered, (b) bond-centered, and (c) asymmetric profiles. 
This progression repeats starting again with the (d) site-centered 
profile until motion ceases due to radiation loss. 

The PN barrier should provide some insight as to where 
these regions may be. While there is no standard formal 
definition of the PN barrier in higher dimensions, one may 
adopt a natural definition (as used in Ref. [2]), according 
to which the barrier is the largest energy difference, for a 
fixed norm of the soliton, between two stationary solutions 
of the system close to configurations that a discrete soli- 
ton must pass when moving adiabatically along the chosen 
lattice direction. This set of configurations includes asym- 
metric solutions, and, importantly in higher dimensions, 
the hybrid solutions too. For example, for a stationary site- 
centered soliton to become mobile along an axis, it must 
overcome barriers created by the asymmetric and hybrid 
states, since, in the course of its motion, its profile will 
change as follows: site-centered — > asymmetric — ► hybrid 
— > asymmetric — > site-centered. If we chose to kick the soli- 
ton along the diagonal, then the same progression should 
be considered with the bond-centered state replacing the 
hybrid one (see Fig. |SJ). Here we use, for the definition of 
the PN barrier, fixed frequency fi rather than fixed norm 
M. We also consider the free energy G = H — fj,M instead 
of the Hamiltonian as in Ref. [29] (note that Eq. (J6j> can be 
derived as dG/dip^.m = 0). 

We kicked the site-centered solutions for various values 
of fc n , and estimated the corresponding threshold values. 
In Fig. [S] the maximum distance traveled, 

Anax(fc)= sup L(n)(t)J - L(n)(0)J, (11) 
te[o,T ] 

where the center of mass is computed by 

(n)(t) =^n|Vv«(*)l7I]l^, m (i)| 2 : (12) 

n.m n,m 

is plotted versus the kicking strength. The corresponding 
threshold values are also identified there. 

It turns out that, the values of the thresholds are related 
to the PN barrier. The left panel in Fig. [TOl shows the dif- 



100r 




Fig. 9. The maximum distance traveled as a function of the kicking 
strength k n for (fj,, C) = (-0.225,0.4) and t G [0,800]. The area 
labeled (a) in the graph represents values of k n that could not depin 
the solution (see Fig. 0a). The area labeled by (b) consists of values 
of k„ that yield a mobile solution (see Fig. 0b) and in (c) the kick 
is so strong that the solution disperses (see Fig. 0c). The threshold 
values, fcjjepjn and fcdisperse are a l so shown. 

ference in free energy, AGhybrid = G s itc — Ghybrid between 
the site-centered solution and the hybrid solution for fixed 
C = 0.4 and \x £ [— 0.3, —0.1]. In each subpanel of the figure 
-Dmax, as defined in Eq. (fTTj) . is plotted against the kicking 
strength for t £ [0, 800] for fixed /U. In panel (i) the site- 
centered solution has more energy than the hybrid solution 
but is unstable and moves away from its initial position even 
for fc„ = 0. Panel (ii) represents parameter values where 
the site-centered solution has greater power and is stable. 
In this small "transparency window" of parameter space, 
there is also pair of unstable asymmetric solutions. In this 
region, we observed the best mobility (see Fig. ITT)). This is 
consistent with what was found in the saturable 2D DNLS 
[14] where good mobility was found where asymmetric so- 
lutions exist. In panel (iii) the threshold fcdisperse is visible 
and the sign of AGhybrid has switched. In (iv) we see that 
the value of fcdepin is increasing and fcdisperse is decreasing 
as the PN barrier increases. Panel (v) corresponds to the 
maximum energy difference. This is also where the largest 
fcdepin occurs. As the energy difference decreases once again 
as seen in (vi) the threshold fcdisperse continues to decrease. 
This is also the case in panel (vii) as both thresholds ap- 
proach fc n = 0. Finally, for the unstable region in (viii) 
fcdepin is once again zero. 

We were unable to identify true transparency points in 
the present model (the 2D DNLS lattice with the CQ on- 
site nonlinearity) , which seems to preclude the possibility 
of finding exact translationally invariant solutions. How- 
ever, enhanced mobility was achieved by lending station- 
ary solutions kinetic energy in cases where the PN barrier 
was low. These moving states gradually lose energy and get 
eventually trapped at some positions in the lattice. Solving 
higher-dimensional counterparts of Eq. ((9]) in the higher- 
dimensional lattice might reveal moving radiationlcss solu- 
tions, although, as we pointed out above, solutions to this 
(quite difficult) problem may not typically exist. It is worth 
mentioning in passing that the energy loss in the ID dis- 



7 




Fig. 10. (Color online) Left: Plot of AGhybrid f° r various values of fj, and fixed G = 0.4. The remaining panels (i)-(viii) correspond to the 
maximum distance traveled versus kicking strength plots. See text for more details. 




AG, 



AG 



AG 




-0.2835 -0.2825 -0.2815 -0.2805 



Fig. 11. (Color online) Top: density plot for the site-centered soliton 
set in motion along the lattice axis for (fi,C) = (—0.282,0.4) and 
k n = 0.5. The choice of parameters fall in a "transparency window" 
where good mobility is observed, possibly due to the existence of a 
pair of asymmetric solutions. A one-dimensional slice along the axis 
of propagation (at m = 10) is shown here. Bottom: zoom of the left 
panel of Fig. 1101 near the "transparency window" . The difference of 
free energy of the site-centered solution and the pair of asymmetric 



solutions AG; 



asymm 



G s 



G a 



is also shown. The energy 



added from the kick exceeds both of these differences. 

Crete sine-Gordon lattice has been recently described using 
an averaged Lagrangian approach in Ref. [40] . 

5. Three-dimensional solutions 

We will now briefly consider a 3D version of the CQ 
DNLS model. The respective counterpart of Eq. ((5]) is 



n,ra,L \ ^r, 



(13) 



where ip n ,m,i is the complex field at site {n,m,l}. In an 
isotropic medium, the discrete Laplacian is 

A (3 Vn,m,; = 1pn+l,m,l + 4>n-l,rn,l + 4>n,m+l,l + 4>n,m~l,l 

+4' ~Hn, m ,l- (14) 



Wc search for stationary solutions, ipn 



rexp(-z^), 



using the same method as in Sec. [2] The 2D soliton species 
have their natural 3D counterparts. As shown in Fig. [T2l 
the extra dimension admits an additional type of a hybrid 
soliton. 





m -2 -2 



Fig. 12. (Color online) Plot of the basic configurations in the 3D 
lattices using iso-contours. Top: plot of 3D site-centered (left) and 
bond-centered (right) solitons. Larger diamonds correspond to larger 
local amplitudes. Bottom: Two different types of 3D hybrid solutions. 
The different colors (arranged in a 3D check-board pattern) are solely 
used for clarity of presentation. 

Figure [TBI shows M(/x) curves for 3D bond-centered and 
site-centered solitons for C = 0.1 and C = 0.7. The figure 
illustrates that in the 3D case, similarly to the 2D case, the 
snake-like patterns are present for small values of coupling 
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Fig. 13. (Color online) The power of the site- and bond-centered 
solitons versus the frequency for C = 0.1 (top) and C = 0.7 (bottom) 
in the 3D lattice. 

constant C and are stretched as C is increased. Similar 
results were obtained for the 3D hybrid solutions (results 
not shown here) . 

6. Variational approximation 

Following the pattern of the VA developed in Ref. [24] 
for ID discrete solitons in the CQ-DNLS model, it is possi- 
ble to construct analytical approximations for the discrete 
solitons, and compare them to the numerical solutions de- 
scribed above. We present this approach for the 2D model, 
but the procedure is essentially the same in three dimen- 
sions. It is relevant to mention that the VA for ID discrete 
solitons in models of the DNLS type was first developed in 
Ref. 05 . 

Solutions to the stationary version Eq. |T]) are local ex- 
trema of the corresponding Lagrangian, 



L ■■ 



E 

n.rn—~ 

C[(t 



1 



(15) 



r n+l,m 



>2 



+ (u. 



n,m+l 



-^n,m) j 



[recall ipn.m, = Un,m exp(— ifit)]. We approximate each soli- 
ton by a localized ansatz which makes it possible to evaluate 
the infinite sums in Eq. (|15p in an explicit form. First, the 
following ansatz is used for the site-centered (sc) solution: 



Fig. 14. Numerical solutions (solid line) and the variational approx- 
imation (triangles) for the site-centered solitons at C = 0.1 in the 
2D lattice model. The approximation based on ansatz (1161 is able to 
capture subfamilies of tall and short narrow solitons, and the branch 
of short wide solitons too. 



,(sc) 



(3 if m = n = 0, 

^g-aCM+M) otherwise 



(16) 



where A, (3, and a are real constants to be found from the 
Euler-Lagrangc equations, 



dL eS dL eS dL c 



0, 



(17) 



OA da dp 
L e g standing for Lagrangian ([15]) evaluated with ansatz 
(TlGl) . In particular, a is treated here as one of the varia- 
tional parameters, in contrast to the ID case, where it was 
expressed in terms of /i and C by means of a relation ob- 
tained from the consideration of the linearized stationary 
equation for decaying "tails" of the soliton [24] , 




a = 2 - ll/C. 



(18) 



We have observed, based on numerous calculations, that 
treating a as a variational parameter yields the same rela- 
tion for a in both the 2D and 3D models. This is consistent 
with solutions in the continuum model where it is known 
that the factor in the exponential tail is independent of the 
dimensional 

Solutions predicted by the VA based on ansatz (Til)]) pro- 
vide for a good fit to the short and tall narrow solutions and 
the first subfamily of wide short solitons of the site-centered 
type, (see Fig.[]3|. At larger values of C, the VA-predictcd 
solutions depart from the numerical ones, which is not sur- 
prising, as the exponential cusp implied by the ansatz is 
not featured by the discrete solitons in the strong-coupling 
(quasi-continuum) model. 

Other solution types can be approximated by appropri- 
ately modified ansatze. In particular, the bond-centered 
(be) soliton is based on a frame built of four points with 



2 In the continuum model the tail decays as r~ 1 / 2 e~ br in the 2D 
case and as r~ 1 e~ br in the 3D case where the factor b is independent 
of the dimension. 




Site-centered 



Fig. 15. Numerical solutions (solid line) and the variational approxi- 
mation for the bond-centered (squares) and hybrid (circles) solitons 
at C = 0.1 in the 2D lattice model. The approximations based on 
the ansatze given in il9t and (1201 ) respectively are able to capture 
subfamilies of tall and short narrow solitons. 

equal amplitudes (see Fig. Ob), whereas the hybrid (hy) 
soliton has just two points in its frame (see Fig. [3Jc). Ac- 
cordingly, the solitons of these types can be modeled by the 
following modifications of ansatz ([To]) : 



,( bc ) 



^4 e -«(|™|+|™|) 

^ e -a(|m-l|+|n|) 

^4 e -«(|™|+|™-i|) 



m, n G {0, 1} 
if m, n < 
if m > 1 , n < 
if m < 0, n > 1 



(19) 



and 



,(hy) 



ie - a (|m-i|+|»-i|) otherwise 

(3 n = 0, me {0,1} 

Ae - a (\m\ + \n\) if m) | n |<0 
Ae -a(|m-l| + |„|) otherwise 



(20) 



Further analysis demonstrates that the modified ansatze 
produce a good approximation for the short and tall narrow 
solutions at small C but not any of the wide families (see 

Fig. us). 

We were also able to predict complicated bifurcations of 
the system by introducing the appropriately chosen asym- 
metric (asym) ansatz: 



(asym) 



01 


n = 0, m = 







fa 


n = 0, m = 


1 




fa 


n = 1, m = 





(21) 


ih 


n = 1, m = 


1 




Ae~ a{ 


|m-fl+|n-CD otherwise 








Fig. 16. Bifurcations featuring the bond-centered, site-centered, and 
asymmetric solutions for C = 0.22. Numerical solutions (lines) and 
its predicted counterparts using the VA based on the ansatz (12 1 l l 
(markers) are in good agreement. The asymmetric VA solution cap- 
tures the main qualitative features of the M(fi) curve (e.g. the dra- 
matic increase of power around fj, m —0.55) but slightly underesti- 
mates the power at the bifurcation points. 

an asymmetric solution. Therefore we have some idea a 
priori what the asymmetric solutions should look like and 
have chosen ansatz (f2Tj) accordingly. For £ = the ansatz 
has the form of a site-centered solution whereas for £ = 
0.5 it will represent a bond-centered solution. All interme- 
diate values of £ represent asymmetric solutions that are 
somewhere between a site-centered and bond-centered so- 
lution. Indeed, the computed value of £ based on the vari- 
ational approximation starts near £ = 0.5 for parameter 
values where the asymmetric solution is almost connected 
to bond-centered solution, and slowly decreases to ( = 
as we alter the parameters until it collides with the site- 
centered solution (see Fig. IT6|). 

Finally we apply the methods to 3D lattice solitons using 
the following site-centered ansatz 



Oc) 

m.n.l 



(22) 



The intention here is to capture the bifurcations where the 
site-centered and bond-centered solutions are connected via 



(3 if m = n = I = 0, 

^g-aCM+M+ltl) otnerwise 

where, for C small enough, it also works well, see Fig. [TTj 

7. Conclusion 

In this work, we have examined the existence, stability, 
and mobility of discrete solitons in 2D and 3D NLS lattices 
with competing (cubic-quintic, CQ) onsite nonlinearities. 
Some properties of the discrete solitons, such as existence 
of the solutions of tall and short types, each narrow and/or 
wide, resemble properties recently found in discrete soli- 
tons in the ID counterpart of this model [23], as well as 
the 2D properties of models such as the one with the sat- 
urable nonlincarity [14] . We have found pitchfork bifurca- 
tions connecting the site-centered and bond-centered soli- 
tons via unstable asymmetric ones, in contrast with the 
ID model, where the connecting asymmetric solutions were 
stable. Another fundamental soliton species that was stud- 
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t-t. A' 



-1 -0.8 -0.6 -0.4 -0.2 

V- 

Fig. 17. Numerical solution (solid line) and the variational approxi- 
mation (triangles) for the site-centered solitons with C = 0.1 using 
the ansatz given in (1 22 I t in the 3D DNLS lattice with cubic-quintic 
nonlincarities. 

icd in this work, viz., hybrid solutions, exists only in the 
higher-dimensional lattice. We have found, in some regions 
of the parameter space, that the site-centered and bond- 
centered solitons were also connected via the hybrid states. 
At small values of the inter-site coupling constant, C, vari- 
ous types of the 2D and 3D stationary discrete solitons are 
well described by the variational approximation (VA) . 

We have also showed that enhanced mobility of 2D dis- 
crete solitons in the CQ lattice can be realized by imparting 
to them kinetic energy exceeding the PN barrier. Never- 
theless, the moving solitons eventually come to a halt, due 
to the radiation loss. In that connection, we were unable 
to find exact transparency points at which translationally 
invariant solutions would be able to exist. However, look- 
ing for carefully crafted radiationless solutions for moving 
solitons in 2D and 3D lattice models remains a challenging 
open problem. It would also be interesting to study the mo- 
bility of the discrete solitons by means of a dynamical ver- 
sion of the the VA (in the ID model with the cubic onsite 
nonlinearity, a dynamical VA was adapted to the analysis of 
collisions between moving discrete solitons in Ref . [Hj , and 
to capture the stationary site-centered and bond-centered 
solutions with a single ansatz in Ref. [42]). 

Getting back to stationary 2D and 3D discrete solitons in 
the cubic-quintic NLS lattice, remaining topics of interest 
are to search for staggered solitons similarly e.g., to the 
work of Ref. [35] for the cubic lattice, as well as lattice 
solitons with intrinsic vorticity. Thus far, discrete lattice 
solitons and vortices were studied in 2D [44] and 3D DNLS 
equation with the cubic nonlinearity 45 . 
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